---
title: "Replication"
output: pdf_document
header-includes:
  - \usepackage{lscape}
---
\begin{document}
```{r, echo=FALSE}
# Install these packages
library(tidyverse)
library(foreign)
library(ggplot2)
library(MASS)
library(haven)
library(sandwich)
library(lmtest)
library(drc)
library(plm)
library("readxl")
library(fixest)
library(stargazer)
library(gt)
```


```{r}
#Loading in the data
data = read_dta('FINALDATASET.dta')
data = data.frame(data)
#data = group_by(data,gname,country)
```

```{r, results="asis"}
#Summary Statistics
data1 = drop_na(data, attack2, lagcompetition2,  lagfactor,  lagcompfactint2, laganocyn, lagdemyn, laggdppc, lagrpc2, lagcivwar, laglogpop)

stargazer(as.data.frame(data1[c('attack2', 'lagcompetition2',  'lagfactor',  'lagcompfactint2', 'laganocyn', 'lagdemyn', 'laggdppc', 'lagrpc2', 'lagcivwar', 'laglogpop')]), type = 'latex', title = 'Summary Statistics')
```

```{r, results="asis"}
# Table 2 -- no ideological sub-setting
model = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data1)

# Now adding clustered standard errors by groupcode
model_clustered = coeftest(model, vcov = vcovCL, cluster = ~groupcode)

stargazer(model_clustered, type = 'latex', title = "Replication Results for All Groups", covariate.labels = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"))
```


```{r}
#---------First Differences for Table 2-----------------

#LAGCOMPETITION2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

#Every covariate set to its mean or median, except lagcompetition2 = mean + 1sd. 
b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE) + sd(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = tibble(exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))

#LAGFACTOR
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE) + sd(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGCOMPFACTINT2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE) + sd(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGANOCYN
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = 0, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = 1, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGDEMYN
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 0, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 1, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGGDPPC
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGRPC2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE)  + sd(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGCIVWAR
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 0, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 1, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGLOGPOP
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE) + sd(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))
```

```{r, results="asis"}
diffs = cbind(diffs, Variable = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"))
colnames(diffs) = c('All Groups', 'Variable')
rownames(diffs) = diffs$Variable
diffs = diffs[1]
stargazer(diffs, type = 'latex', summary = FALSE, title = 'First Differences for All Groups')
diffs_later = diffs
first_diffs = diffs
```



```{r}
# Table 3, Model 1 -- religious groups
data_religious = data[data$RANDrel==1, ]

model1 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_religious, control=glm.control(maxit=1000), init.theta = 0.1) 

model1_clustered = coeftest(model1, vcov = vcovCL, cluster = ~groupcode)


```
```{r}
#Checking the proportion of zeros in the religious data
#sum(data_religious$attack2 == 0, na.rm=TRUE)/length(data_religious$attack2)
```


```{r}

# Table 3, Model 2 -- nationalist groups
data_nationalist = data[data$RANDnat==1, ]

model2 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_nationalist, control=glm.control(maxit=1000)) 

model2_clustered = coeftest(model2, vcov = vcovCL, cluster = ~groupcode)

```

```{r}
# Table 3, Model 3 -- left-wing groups
data_left = data[data$RANDlw==1, ]

model3 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_left, control=glm.control(maxit=1000)) 

model3_clustered = coeftest(model3, vcov = vcovCL, cluster = ~groupcode)
```

```{r}
# Table 3, Model 4 -- right-wing groups
data_right = data[data$RANDrw==1, ]

model4 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_right, control=glm.control(maxit=1000)) 

model4_clustered = coeftest(model4, vcov = vcovCL, cluster = ~groupcode)
```

```{r, results="hide"}
#Table with the results
stargazer(model1_clustered, model2_clustered, model3_clustered, model4_clustered, type = 'latex', title = "Replication Results by Group Ideology", dep.var.caption = '' , covariate.labels = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), model.numbers = FALSE, column.labels = c('Religious', 'Nationalist', 'Left-Wing ', 'Right-Wing'))
```
\begin{landscape}
\begin{table}[!htbp] \centering 
  \caption{Replication Results by Group Ideology} 
  \label{} 
\begin{tabular}{@{\extracolsep{5pt}}lcccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
\\[-1.8ex] & \multicolumn{4}{c}{ } \\ 
 & Religious & Nationalist & Left-Wing  & Right-Wing \\ 
\hline \\[-1.8ex] 
 Competition & $-$1.261$^{***}$ & $-$0.733$^{***}$ & $-$0.406 & $-$0.321 \\ 
  & (0.439) & (0.168) & (0.257) & (0.370) \\ 
  & & & & \\ 
 Acceptability of Violence & $-$1.525$^{***}$ & $-$0.253 & 0.769$^{**}$ & 1.200$^{**}$ \\ 
  & (0.473) & (0.430) & (0.366) & (0.551) \\ 
  & & & & \\ 
 Competition x Acceptability & 0.725$^{***}$ & 0.094 & $-$0.084 & $-$0.429 \\ 
  & (0.243) & (0.126) & (0.172) & (0.404) \\ 
  & & & & \\ 
 Anocracy & 0.153 & 2.440$^{***}$ & $-$0.711 & $-$2.982$^{***}$ \\ 
  & (1.094) & (0.478) & (0.455) & (0.746) \\ 
  & & & & \\ 
 Democracy & 0.497 & 2.442$^{***}$ & 1.113$^{*}$ & 0.096 \\ 
  & (1.684) & (0.429) & (0.568) & (0.318) \\ 
  & & & & \\ 
 GDP & $-$0.0001$^{*}$ & 0.0001$^{***}$ & 0.00001 & 0.00001 \\ 
  & (0.00005) & (0.00003) & (0.00003) & (0.00005) \\ 
  & & & & \\ 
 RPC & $-$0.462 & $-$0.484 & $-$0.797 & 0.486$^{*}$ \\ 
  & (0.905) & (0.477) & (0.517) & (0.276) \\ 
  & & & & \\ 
 Civil War & 1.323 & 0.716$^{**}$ & 0.474 & 0.949$^{**}$ \\ 
  & (1.379) & (0.327) & (0.416) & (0.386) \\ 
  & & & & \\ 
 Logged Population & $-$0.260 & 0.290$^{**}$ & $-$0.631$^{***}$ & $-$0.391$^{*}$ \\ 
  & (0.370) & (0.118) & (0.117) & (0.230) \\ 
  & & & & \\ 
 Constant & 4.984 & $-$2.673 & 8.483$^{***}$ & 3.735 \\ 
  & (4.424) & (1.796) & (1.323) & (2.375) \\ 
  & & & & \\ 
\hline \\[-1.8ex] 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{4}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\ 
\end{tabular} 
\end{table}  
\end{landscape}

```{r}
#---------First Differences for Table 3-----------------

models = list(model1, model2, model3, model4)

for(model in models) {

#LAGCOMPETITION2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

#Every covariate set to its mean or median, except lagcompetition2 = mean + 1sd. 
b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE) + sd(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = tibble(exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))

#LAGFACTOR
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE) + sd(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGCOMPFACTINT2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE) + sd(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGANOCYN
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = 0, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = 1, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGDEMYN
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 0, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 1, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGGDPPC
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))


#LAGRPC2
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE)  + sd(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGCIVWAR
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 0, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 1, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))



#LAGLOGPOP
a =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('lagcompetition2' = mean(data$lagcompetition2, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactint2' = mean(data$lagcompfactint2, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE) + sd(data$laglogpop, na.rm = TRUE))

diffs = rbind(diffs, exp(predict(model, newdata=b)) - exp(predict(model, newdata=a)))
first_diffs = cbind(first_diffs, diffs)
}
colnames(first_diffs) = c('All Groups', 'Religious', 'Nationalist', 'Left-Wing', 'Right-Wing')
```

```{r, results="hide"}
stargazer(first_diffs, type = 'latex', summary = FALSE, title = "First Differences by Group Ideology", dep.var.caption = '')
```
\begin{landscape}
\begin{table}[!htbp] \centering 
  \caption{First Differences by Group Ideology} 
  \label{} 
\begin{tabular}{@{\extracolsep{5pt}} cccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
 & All Groups & Religious & Nationalist & Left-Wing & Right-Wing \\ 
\hline \\[-1.8ex] 
Competition & $$-$7.055$ & $$-$1.969$ & $$-$4.639$ & $$-$4.253$ & $$-$0.796$ \\ 
Acceptability of Violence & $7.558$ & $$-$2.099$ & $$-$1.981$ & $15.119$ & $7.035$ \\ 
Competition x Acceptability & $0.322$ & $12.268$ & $2.087$ & $$-$2.128$ & $$-$1.690$ \\ 
Anocracy & $2.478$ & $0.430$ & $86.316$ & $$-$5.892$ & $$-$2.491$ \\ 
Democracy & $10.357$ & $1.015$ & $7.528$ & $7.774$ & $0.239$ \\ 
GDP & $5.904$ & $$-$0.973$ & $4.775$ & $0.535$ & $0.215$ \\ 
RPC & $$-$2.621$ & $$-$0.561$ & $$-$1.858$ & $$-$3.972$ & $0.765$ \\ 
Civil War & $12.404$ & $0.930$ & $8.992$ & $4.909$ & $1.824$ \\ 
Logged Population & $$-$5.664$ & $$-$0.873$ & $4.778$ & $$-$7.296$ & $$-$1.206$ \\ 
\hline \\[-1.8ex] 
\end{tabular} 
\end{table} 
\end{landscape}

```{r, results="hide"}
#Table with the results
stargazer(model1_clustered, model2_clustered, model3_clustered, model4_clustered, type = 'latex', title = "Replication Results by Group Ideology", dep.var.caption = '' , covariate.labels = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), model.numbers = FALSE, column.labels = c('Religious', 'Nationalist', 'Left-Wing ', 'Right-Wing'))

```

# Intervention 1 Table 2, a new competition variable // competition_delta for lagcompetition2
```{r, results='asis'}
# Intervention 1 Table 2, a new competition variable // competition_delta for lagcompetition2
data_dplyr <- data %>%                 
  group_by(groupcode,country) %>%
  dplyr::mutate(lag2 = dplyr::lag(lagcompetition2, n = 1, default = NA)) %>% 
  as.data.frame()

data <- data_dplyr %>%                            
  group_by(groupcode,country) %>%
  dplyr::mutate(lag3 = dplyr::lag(lag2, n = 1, default = NA)) %>% 
  as.data.frame()

data$competition_delta = data$lagcompetition2 - data$lag2
data$lagcompfactor_delta = data$competition_delta * data$lagfactor

model_int = glm.nb(attack2 ~ competition_delta + lagfactor  + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data, control=glm.control(maxit=1000))

# Now adding clustered standard errors by groupcode
model_int_clustered = coeftest(model_int, vcov = vcovCL, cluster = ~groupcode)

#results table
stargazer(model_clustered, model_int_clustered, type = 'text', title = "Intervention 1, All Groups", covariate.labels = c('Nameth Competition', 'Intervention Revised Competition', "Acceptability of Violence", "Nameth x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), dep.var.caption = '')

```
```{r}
#---------First Differences for I1, Table 2-----------------



#COMPETITION_DELTA
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

#Every covariate set to its mean or median, except competition_delta = mean + 1sd. 
b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE) + sd(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = tibble(exp(predict(model_int, newdata=b) - exp(predict(model_int, newdata=a))))


#LAGFACTOR
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE) + sd(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs_int, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))


#LAGCOMPFACTOR_DELTA
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE) + sd(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))


#LAGANOCYN
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = 0, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = 1, 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))


#LAGDEMYN
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 0, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = 1, 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))



#LAGGDPPC
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))


#LAGRPC2
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE)  + sd(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))



#LAGCIVWAR
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 0, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE) + sd(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = 1, 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))



#LAGLOGPOP
a =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE))

b =  data.frame('competition_delta' = mean(data$competition_delta, na.rm = TRUE), 'lagfactor' = mean(data$lagfactor, na.rm = TRUE), 'lagcompfactor_delta' = mean(data$lagcompfactor_delta, na.rm = TRUE), 'laganocyn' = median(data$laganocyn, na.rm = TRUE), 'lagdemyn' = median(data$lagdemyn, na.rm = TRUE), 'laggdppc' = mean(data$laggdppc, na.rm = TRUE), 'lagrpc2' = mean(data$lagrpc2, na.rm = TRUE), 'lagcivwar' = median(data$lagcivwar, na.rm = TRUE), 'laglogpop' = mean(data$laglogpop, na.rm = TRUE) + sd(data$laglogpop, na.rm = TRUE))

diffs_int = rbind(diffs, exp(predict(model_int, newdata=b)) - exp(predict(model_int, newdata=a)))
```

```{r, results="asis"}
diffs_int = cbind(diffs, Variable = c('Competition ', "Acceptability of Violence", "Intervention x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"))
colnames(diffs_int) = c('All Groups', 'Variable')
rownames(diffs_int) = diffs_int$Variable
diffs_int = diffs_int[1]

diffs_int = cbind(diffs_later, diffs_int)
colnames(diffs_int) = c('Nameth Competition', 'Intervention Competition')
#stargazer(diffs_int, type = 'latex', summary = FALSE, title = 'Competition Measure Intervention, All Groups')

```


```{r}
# I1 Table 3, Model 1
data_religious = data[data$RANDrel==1, ]

model1_int = glm.nb(attack2 ~ competition_delta + lagfactor + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_religious, control=glm.control(maxit=1000)) 

model1_clustered_int = coeftest(model1_int, vcov = vcovCL, cluster = ~groupcode)
```
```{r}
# I1 Table 3, Model 2 -- nationalist groups
data_nationalist = data[data$RANDnat==1, ]

model2_int = glm.nb(attack2 ~ competition_delta + lagfactor + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_nationalist, control=glm.control(maxit=1000)) 

model2_clustered_int = coeftest(model2_int, vcov = vcovCL, cluster = ~groupcode)
```
```{r}
# I1 Table 3, Model 3 -- left-wing groups
data_left = data[data$RANDlw==1, ]

model3_int = glm.nb(attack2 ~ competition_delta + lagfactor + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_left, control=glm.control(maxit=1000)) 

model3_clustered_int = coeftest(model3_int, vcov = vcovCL, cluster = ~groupcode)
```
```{r}
# I1 Table 3, Model 4 -- right-wing groups
data_right = data[data$RANDrw==1, ]

model4_int = glm.nb(attack2 ~ competition_delta + lagfactor  + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_right, control=glm.control(maxit=1000)) 

model4_clustered_int = coeftest(model4_int, vcov = vcovCL, cluster = ~groupcode)
```

```{r, results = 'hide'}
#Results

stargazer(model1_clustered, model1_clustered_int, model2_clustered, model2_clustered_int, model3_clustered, model3_clustered_int, model4_clustered, model4_clustered_int, type = 'latex', title = "Competition Measure Intervention by Group Ideology", dep.var.caption = '' , model.numbers = FALSE, column.labels = c('Religious (N)','Religious (I)' , 'Nationalist (N)', 'Nationalist (I)', 'Left-Wing (N)', 'Left-Wing (I)', 'Right-Wing (N)', 'Right-Wing (I)'), covariate.labels = c('Nameth Competition', 'Intervention Revised Competition', "Acceptability of Violence", "Nameth x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), notes = 'The original models are represented with an I, our results are represented with an N', notes.align = 'r', font.size = 'tiny')
```
\begin{landscape}
\begin{table}[!htbp] \centering 
  \caption{Intervention 1, by Ideology} 
  \label{} 
\tiny 
\begin{tabular}{@{\extracolsep{5pt}}lcccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
\\[-1.8ex] & \multicolumn{8}{c}{ } \\ 
 & Religious (N) & Religious (I) & Nationalist (N) & Nationalist (I) & Left-Wing (N) & Left-Wing (I) & Right-Wing (N) & Right-Wing (I) \\ 
\hline \\[-1.8ex] 
 Nameth Competition & $-$1.261$^{***}$ &  & $-$0.733$^{***}$ &  & $-$0.406 &  & $-$0.321 &  \\ 
  & (0.439) &  & (0.168) &  & (0.257) &  & (0.370) &  \\ 
  & & & & & & & & \\ 
 Intervention Revised Competition &  & $-$0.204$^{***}$ &  & 0.019 &  & 0.222$^{***}$ &  & 0.102 \\ 
  &  & (0.057) &  & (0.025) &  & (0.052) &  & (0.231) \\ 
  & & & & & & & & \\ 
 Acceptability of Violence & $-$1.525$^{***}$ & $-$1.028$^{***}$ & $-$0.253 & $-$0.307 & 0.769$^{**}$ & 0.492$^{***}$ & 1.200$^{**}$ & 0.720$^{***}$ \\ 
  & (0.473) & (0.337) & (0.430) & (0.289) & (0.366) & (0.174) & (0.551) & (0.267) \\ 
  & & & & & & & & \\ 
 Nameth x Acceptability & 0.725$^{***}$ &  & 0.094 &  & $-$0.084 &  & $-$0.429 &  \\ 
  & (0.243) &  & (0.126) &  & (0.172) &  & (0.404) &  \\ 
  & & & & & & & & \\ 
 Anocracy & 0.153 & $-$0.160 & 2.440$^{***}$ & 2.340$^{***}$ & $-$0.711 & $-$0.670$^{*}$ & $-$2.982$^{***}$ & $-$2.897$^{***}$ \\ 
  & (1.094) & (1.183) & (0.478) & (0.625) & (0.455) & (0.403) & (0.746) & (0.747) \\ 
  & & & & & & & & \\ 
 Democracy & 0.497 & 0.315 & 2.442$^{***}$ & 2.243$^{***}$ & 1.113$^{*}$ & 0.828 & 0.096 & $-$0.077 \\ 
  & (1.684) & (1.681) & (0.429) & (0.492) & (0.568) & (0.526) & (0.318) & (0.358) \\ 
  & & & & & & & & \\ 
 GDP & $-$0.0001$^{*}$ & $-$0.0001$^{**}$ & 0.0001$^{***}$ & 0.0001$^{**}$ & 0.00001 & $-$0.00000 & 0.00001 & 0.00001 \\ 
  & (0.00005) & (0.00005) & (0.00003) & (0.00003) & (0.00003) & (0.00003) & (0.00005) & (0.0001) \\ 
  & & & & & & & & \\ 
 RPC & $-$0.462 & $-$0.514 & $-$0.484 & $-$0.538 & $-$0.797 & $-$0.753 & 0.486$^{*}$ & 0.432 \\ 
  & (0.905) & (1.023) & (0.477) & (0.518) & (0.517) & (0.509) & (0.276) & (0.288) \\ 
  & & & & & & & & \\ 
 Civil War & 1.323 & 1.349 & 0.716$^{**}$ & 0.755$^{*}$ & 0.474 & 0.660 & 0.949$^{**}$ & 1.185$^{***}$ \\ 
  & (1.379) & (1.303) & (0.327) & (0.416) & (0.416) & (0.424) & (0.386) & (0.409) \\ 
  & & & & & & & & \\ 
 Logged Population & $-$0.260 & $-$0.409 & 0.290$^{**}$ & 0.512$^{***}$ & $-$0.631$^{***}$ & $-$0.711$^{***}$ & $-$0.391$^{*}$ & $-$0.297 \\ 
  & (0.370) & (0.372) & (0.118) & (0.139) & (0.117) & (0.126) & (0.230) & (0.331) \\ 
  & & & & & & & & \\ 
 Constant & 4.984 & 5.827 & $-$2.673 & $-$5.692$^{***}$ & 8.483$^{***}$ & 8.915$^{***}$ & 3.735 & 2.495 \\ 
  & (4.424) & (4.801) & (1.796) & (2.105) & (1.323) & (1.427) & (2.375) & (3.312) \\ 
  & & & & & & & & \\ 
\hline \\[-1.8ex] 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{8}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\ 
 & \multicolumn{8}{r}{The original models are represented with an I, our results are represented with an N} \\ 
\end{tabular} 
\end{table} 
\end{landscape}




# Intervention 2, introducing country fixed effects
```{r, results='asis'}
# Intervention 2, introducing country fixed effects
model = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data, control=glm.control(maxit=1000)) 

fe.model <- feglm(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop | factor(country), data = data, family = negative.binomial(model$theta)) ## this is the dispersion parameter provided by summary(model)

# Now adding clustered standard errors by groupcode AND calculating variance-covariance matrix for feglm model.
fe.vcovCL = vcov(fe.model)
model_clustered = coeftest(fe.model, vcov = fe.vcovCL, cluster = ~groupcode)

stargazer(model, model_clustered, type = 'latex', title = "Fixed Effects Intervention, All Groups", dep.var.labels.include = FALSE , model.numbers = FALSE, covariate.labels = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), model.names = FALSE, column.labels = c('Original Model', 'FE Model'), dep.var.caption = "")
```

```{r}
# I3 Table 3, Model 1 -- religious groups
data_religious = data[data$RANDrel==1, ]

model1 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_religious, control=glm.control(maxit=1000)) 

model1_clustered = coeftest(model1, vcov = vcovCL, cluster = ~groupcode)

fe.model1 <- feglm(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop | factor(country), data = data_religious, family = negative.binomial(model1$theta)) ## this is the dispersion parameter provided by summary(model1)

# Now adding clustered standard errors by groupcode AND calculating variance-covariance matrix for feglm model.
fe.vcovCL1 = vcov(fe.model1)
model1_clustered_int = coeftest(fe.model1, vcov = fe.vcovCL1, cluster = ~groupcode)
```

```{r}
# I3 Table 3, Model 2 -- nationalist groups
data_nationalist = data[data$RANDnat==1, ]

model2 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_nationalist, control=glm.control(maxit=1000)) 

model2_clustered = coeftest(model2, vcov = vcovCL, cluster = ~groupcode)

fe.model2 <- feglm(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop | factor(country), data = data_nationalist, family = negative.binomial(model2$theta)) ## this is the dispersion parameter provided by summary(model2)

# Now adding clustered standard errors by groupcode AND calculating variance-covariance matrix for feglm model.
fe.vcovCL2 = vcov(fe.model2)
model2_clustered_int = coeftest(fe.model2, vcov = fe.vcovCL2, cluster = ~groupcode)
```


```{r}
# I3 Table 3, Model 3 -- left-wing groups
data_left = data[data$RANDlw==1, ]

model3 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_left, control=glm.control(maxit=1000)) 

model3_clustered = coeftest(model3, vcov = vcovCL, cluster = ~groupcode)

fe.model3 <- feglm(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop | factor(country), data = data_left, family = negative.binomial(model3$theta)) ## this is the dispersion parameter provided by summary(model2)

# Now adding clustered standard errors by groupcode AND calculating variance-covariance matrix for feglm model.
fe.vcovCL3 = vcov(fe.model3)
model3_clustered_int = coeftest(fe.model3, vcov = fe.vcovCL3, cluster = ~groupcode)
```

```{r}
# I3 Table 3, Model 4 -- right-wing groups
data_right = data[data$RANDrw==1, ]

model4 = glm.nb(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop, data = data_right, control=glm.control(maxit=1000))

model4_clustered = coeftest(model4, vcov = vcovCL, cluster = ~groupcode)

fe.model4 <- feglm(attack2 ~ lagcompetition2 + lagfactor + lagcompfactint2 + laganocyn + lagdemyn + laggdppc + lagrpc2 + lagcivwar + laglogpop | factor(country), data = data_right, family = negative.binomial(model3$theta)) ## this is the dispersion parameter provided by summary(model2)

# Now adding clustered standard errors by groupcode AND calculating variance-covariance matrix for feglm model.
fe.vcovCL4 = vcov(fe.model4)
model4_clustered_int = coeftest(fe.model4, vcov = fe.vcovCL4, cluster = ~groupcode)
```

```{r, results = 'hide'}
#Table with the results
stargazer(model1_clustered, model1_clustered_int, model2_clustered, model2_clustered_int, model3_clustered, model3_clustered_int, model4_clustered, model4_clustered_int, type = 'latex', title = "Fixed Effects Intervention by Group Ideology", dep.var.caption = '' , covariate.labels = c('Competition', "Acceptability of Violence", "Competition x Acceptability", "Anocracy", "Democracy", "GDP", "RPC", "Civil War", "Logged Population"), model.numbers = FALSE, column.labels = c('Religious (N)','Religious (I)' , 'Nationalist (N)', 'Nationalist (I)', 'Left-Wing (N)', 'Left-Wing (I)', 'Right-Wing (N)', 'Right-Wing (I)'), font.size='tiny')
```
\begin{landscape}
\begin{table}[!htbp] \centering 
  \caption{Fixed Effects Intervention by Group Ideology} 
  \label{} 
\tiny 
\begin{tabular}{@{\extracolsep{5pt}}lcccccccc} 
\\[-1.8ex]\hline 
\hline \\[-1.8ex] 
\\[-1.8ex] & \multicolumn{8}{c}{ } \\ 
 & Religious (N) & Religious (I) & Nationalist (N) & Nationalist (I) & Left-Wing (N) & Left-Wing (I) & Right-Wing (N) & Right-Wing (I) \\ 
\hline \\[-1.8ex] 
 Competition & $-$1.261$^{***}$ & $-$0.911 & $-$0.733$^{***}$ & $-$0.230$^{***}$ & $-$0.406 & $-$0.300 & $-$0.321 & $-$0.089 \\ 
  & (0.439) & (1.029) & (0.168) & (0.078) & (0.257) & (0.187) & (0.370) & (0.425) \\ 
  & & & & & & & & \\ 
 Acceptability of Violence & $-$1.525$^{***}$ & $-$0.640 & $-$0.253 & 0.194 & 0.769$^{**}$ & 0.453 & 1.200$^{**}$ & 0.064 \\ 
  & (0.473) & (0.904) & (0.430) & (0.197) & (0.366) & (0.289) & (0.551) & (0.654) \\ 
  & & & & & & & & \\ 
 Competition x Acceptability & 0.725$^{***}$ & 0.696 & 0.094 & $-$0.012 & $-$0.084 & $-$0.084 & $-$0.429 & $-$0.165 \\ 
  & (0.243) & (0.608) & (0.126) & (0.065) & (0.172) & (0.151) & (0.404) & (0.519) \\ 
  & & & & & & & & \\ 
 Anocracy & 0.153 & $-$1.298$^{***}$ & 2.440$^{***}$ & 1.166$^{***}$ & $-$0.711 & $-$0.292 & $-$2.982$^{***}$ & $-$2.917$^{***}$ \\ 
  & (1.094) & (0.428) & (0.478) & (0.365) & (0.455) & (0.371) & (0.746) & (0.503) \\ 
  & & & & & & & & \\ 
 Democracy & 0.497 & $-$1.554$^{***}$ & 2.442$^{***}$ & 2.021$^{***}$ & 1.113$^{*}$ & 0.879$^{*}$ & 0.096 & 0.544 \\ 
  & (1.684) & (0.578) & (0.429) & (0.391) & (0.568) & (0.475) & (0.318) & (0.684) \\ 
  & & & & & & & & \\ 
 GDP & $-$0.0001$^{*}$ & $-$0.0002$^{*}$ & 0.0001$^{***}$ & $-$0.0001$^{**}$ & 0.00001 & $-$0.0001$^{**}$ & 0.00001 & $-$0.0001 \\ 
  & (0.00005) & (0.0001) & (0.00003) & (0.00002) & (0.00003) & (0.0001) & (0.00005) & (0.0001) \\ 
  & & & & & & & & \\ 
 RPC & $-$0.462 & $-$1.259 & $-$0.484 & $-$0.338 & $-$0.797 & $-$1.087$^{**}$ & 0.486$^{*}$ & 0.628$^{***}$ \\ 
  & (0.905) & (0.773) & (0.477) & (0.339) & (0.517) & (0.480) & (0.276) & (0.036) \\ 
  & & & & & & & & \\ 
 Civil War & 1.323 & 1.353$^{*}$ & 0.716$^{**}$ & 0.991$^{***}$ & 0.474 & 0.586 & 0.949$^{**}$ & 1.532$^{*}$ \\ 
  & (1.379) & (0.755) & (0.327) & (0.312) & (0.416) & (0.361) & (0.386) & (0.874) \\ 
  & & & & & & & & \\ 
 Logged Population & $-$0.260 & 7.818$^{***}$ & 0.290$^{**}$ & 2.829$^{***}$ & $-$0.631$^{***}$ & 1.056 & $-$0.391$^{*}$ & 4.293$^{*}$ \\ 
  & (0.370) & (2.312) & (0.118) & (1.000) & (0.117) & (0.946) & (0.230) & (2.593) \\ 
  & & & & & & & & \\ 
 Constant & 4.984 &  & $-$2.673 &  & 8.483$^{***}$ &  & 3.735 &  \\ 
  & (4.424) &  & (1.796) &  & (1.323) &  & (2.375) &  \\ 
  & & & & & & & & \\ 
\hline \\[-1.8ex] 
\hline 
\hline \\[-1.8ex] 
\textit{Note:}  & \multicolumn{8}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\ 
\end{tabular} 
\end{table} 
\end{landscape}



# Figures
```{r}
#Figure 1, observations per country
no_obs = data %>% 
  group_by(country) %>% 
  summarise(n = n())
figure1 = ggplot(no_obs, aes(x=reorder(country, -n), y=n)) + geom_bar(stat="identity") +labs(title="Number of Observations by Country",
        x ="Countries", y = "No. of Observations") + theme(plot.title = element_text(hjust = 0.5), axis.text.x=element_blank(),panel.background = element_blank(), axis.ticks.x=element_blank(),panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank())
figure1

#Figure 2, average attacks per year per country
avg_attack = data %>%
  group_by(country,year) %>%
  summarise_at(vars(attack2), list(sum = sum), na.rm=TRUE)
avg_attack = avg_attack %>%
  group_by(country) %>%
  summarise_at(vars(sum), list(avg = mean), na.rm=TRUE)
avg_attack

figure2 = ggplot(avg_attack, aes(x=reorder(country, -avg), y=avg)) + geom_bar(stat="identity") +labs(title="Average Annual Attacks by Country",
        x ="Countries", y = "Mean of Group Attacks") + theme(plot.title = element_text(hjust = 0.5), axis.text.x=element_blank(),panel.background = element_blank(), axis.ticks.x=element_blank(),panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank())
figure2
#write.csv(avg_attack, file="avg_attack.csv")

#Figure 3, average competition index per year per country
avg_comp = data %>%
  group_by(country) %>%
  summarise_at(vars(lagcompetition2), list(avg = mean), na.rm=TRUE)
avg_comp
avg_comp = data %>%
  group_by(country) %>%
  summarise_at(vars(lagcompetition2), list(avg = mean), na.rm=TRUE)
avg_comp
#write.csv(avg_comp, file="avg_comp.csv")

figure3 = ggplot(avg_comp, aes(x=reorder(country, -avg), y=avg)) + geom_bar(stat="identity") +labs(title="Average Annual Competition Index by Country",
        x ="Countries", y = "Mean Competition Index Value") + theme(plot.title = element_text(hjust = 0.5), axis.text.x=element_blank(),panel.background = element_blank(), axis.ticks.x=element_blank(),panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank())
figure3

#Figure 3, histogram of annual competition index by country
comp_hist = data[!is.na(data$lagcompetition2),]
comp_hist = data %>% distinct(country, year, .keep_all = TRUE)

figure4 = ggplot(comp_hist, aes(x=lagcompetition2)) + geom_histogram() +stat_bin(binwidth = 1) +labs(title="Country-Year Competition Index Values",
        x ="Competition Index Values", y = "Observations") + theme(plot.title = element_text(hjust = 0.5), panel.background = element_blank(), axis.ticks.x=element_blank(),panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank())
figure4

```



\end{document}
